
typedef struct 
{
    float4 pos;
    float4 color;
    float4 particleProp;
} ParticleRender;

typedef struct
{
    float4 speed;
} ParticleSim;

const sampler_t samplerLinear = CLK_FILTER_LINEAR | CLK_NORMALIZED_COORDS_TRUE; 
const sampler_t samplerLinear2 = CLK_FILTER_LINEAR | CLK_NORMALIZED_COORDS_FALSE; 
const sampler_t samplerNearest = CLK_FILTER_NEAREST | CLK_NORMALIZED_COORDS_TRUE; 

kernel void firstKernel(    global ParticleRender* renderPtr, global ParticleSim* simPtr, 
                            read_only image2d_t input, write_only image2d_t output,
                            read_only image2d_t normalTex,
                            unsigned int offset, unsigned int count,
                            float dt,
                            unsigned int texWidth, unsigned int texHeight,
                            float terrainWidth, float terrainHeight)
{
    //const float3 g = (float3)(0.0, 0.0, 0.0);
    const float3 g = (float3)(0.0, -0.0075, 0.0);
    const float damp = 0.98;

    unsigned int idx = get_global_id(0);

    if ( idx >= offset+count || renderPtr[idx].particleProp.y >= renderPtr[idx].particleProp.z)
        return;

    /* get particle properties */
    float4 pos  = renderPtr[idx].pos;
    float3 v    = simPtr[idx].speed.xyz;
        
    /* process */ 

    // integrate 
    v += g;
    pos.xyz += v;

    //float2 tex2 = (float2)(texWidth*pos.x/terrainWidth, texHeight*pos.z/terrainHeight); //[0,texDim]

    float2 tex = (float2)(pos.x/terrainWidth, pos.z/terrainHeight); // [0,1]
    //float2 normTex = (float2)(tex.x/texWidth, tex.y/texHeight);
    //int2 iTex = (int2)(floor(tex.x*(texWidth-1)), floor(tex.y*(texHeight-1)));
    //float h = read_imagef(input, samplerLinear, tex).x;
    //float h0 = read_imagef(input, samplerLinear, (float2)((iTex.x+0.5)/texWidth, (iTex.y+0.5)/texHeight)).x;
    float h0 = read_imagef(input, samplerLinear, tex + (float2)(0.5/texWidth, 0.5/texHeight)).x;

    float3 normal = -read_imagef(normalTex, samplerLinear, tex + (float2)(0.5/texWidth, 0.5/texHeight)).xyz;
    // are we really close to the terrain 
    float radius = renderPtr[idx].particleProp.x;
    if ( (pos.y - h0) < (radius + 1e-5) ) 
    {
        pos.y = radius + h0 + 1e-4; // adjust so were slightly over the terrain 
        //v.y *= -1.0 * damp; 
        v *= damp;
        float3 vN = normalize(v);
        float3 refl = vN - 2.0*dot(vN, normal)*normal;
        v = normalize(refl) * length(v);
        //v = normal * length(v);
#if 1
        float2 fTex = (float2)(tex.x*texWidth, tex.y*texHeight); // [0.0, texDim-1]
        int2 iTex = (int2)(tex.x*texWidth, tex.y*texHeight);     // [0,texDim-1]

        float offsetX = fTex.x - iTex.x;
        float offsetY = fTex.y - iTex.y;

#if 0
        float d0 = sqrt((1.0-offsetX)*(1.0-offsetX) + (1.0-offsetY)*(1.0-offsetY));
        float d1 = sqrt(offsetX*offsetX             + (1.0-offsetY)*(1.0-offsetY));
        float d2 = sqrt((1.0-offsetX)*(1.0-offsetX) + offsetY*offsetY);
        float d3 = sqrt(offsetX*offsetX             + offsetY*offsetY);
#else
        float d0 = 1.41421356 - sqrt(offsetX*offsetX             + offsetY*offsetY);
        float d1 = 1.41421356 - sqrt((1.0-offsetX)*(1.0-offsetX) + offsetY*offsetY);
        float d2 = 1.41421356 - sqrt(offsetX*offsetX             + (1.0-offsetY)*(1.0-offsetY));
        float d3 = 1.41421356 - sqrt((1.0-offsetX)*(1.0-offsetX) + (1.0-offsetY)*(1.0-offsetY));
        //float d0 = (1.0-offsetX) + (1.0-offsetY);
        //float d1 = offsetX       + (1.0-offsetY);
        //float d2 = (1.0-offsetX) + offsetY);
        //float d3 = offsetX       + offsetY;
#endif
        float hn0 = read_imagef(input, samplerNearest, tex + (float2)(0.0, 0.0)).x;
        float hn1 = read_imagef(input, samplerNearest, tex + (float2)(1.0/texWidth, 0.0)).x;
        float hn2 = read_imagef(input, samplerNearest, tex + (float2)(0.0, 1.0/texHeight)).x;
        float hn3 = read_imagef(input, samplerNearest, tex + (float2)(1.0/texWidth, 1.0/texHeight)).x;

        float s = 0.05;
#if 1 // enabled terrain deformation 
#if 1
        if ( iTex.x > 0 && iTex.y > 0 && iTex.x < texWidth-1 && iTex.y < texHeight-1)
        {
            write_imagef(output, iTex + (int2)(0,0), hn0-d0*s);    
            write_imagef(output, iTex + (int2)(1,0), hn1-d1*s);
            write_imagef(output, iTex + (int2)(0,1), hn2-d2*s);
            write_imagef(output, iTex + (int2)(1,1), hn3-d3*s);
        }
#else
        
        write_imagef(output, iTex, h0-s);   
#endif
#endif
#endif 
    }

    float4 color = renderPtr[idx].color;

    /* write back */     
    simPtr[idx].speed     = (float4)(v.xyz,0.0);
    renderPtr[idx].pos    = pos;

    float4 col0 = (float4)(0.0, 1.0, 0.0, 1.0);
    float4 col1 = (float4)(1.0, 0.0, 0.0, 1.0);

    float alpha = min( fabs(v.y*3.0), 1.0);

    //if ( v.y > 0.0)
    //    renderPtr[idx].color = col0;
    //else
    //    renderPtr[idx].color = col1;

    renderPtr[idx].particleProp.y += dt;
    //renderPtr[idx].color = (1.0-alpha)*col0 + alpha*col1;
    //renderPtr[idx].color = smoothstep(col0, col1, 1.0);
}

kernel void copyKernel(     read_only image2d_t input, write_only image2d_t output,
                            unsigned int width, unsigned int height)
{
    int2 pos = (int2)(get_global_id(0), get_global_id(1));
    float2 tex = (float2)((float)pos.x/width, (float)pos.y/height);

    if(pos.x >= width || pos.y >= height)
        return;

    float4 value = read_imagef(input, samplerNearest, tex);
    write_imagef(output, pos, value);
}
